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An equilibrium random surface multistep height model proposed in [Abraham and Newman, EPL, 

86, 16002 (2009)] is studied using a variant of the worm algorithm. In one limit, the model reduces 

to the two-dimensional Ising model in the height representation. When the Ising model constraint 

of single height steps is relaxed, the critical temperature and critical exponents are continuously 

£_) varying functions of the parameter controlling height steps larger than one. Numerical estimates 

of the critical exponents can be mapped via a single parameter- the Coulomb gas coupling- to the 

exponents of the O(n) loop model on the honeycomb lattice with n < 1. 
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CO 
+J I. INTRODUCTION 

d 

Some time ago, two of the authors introduced a statistical mechanical model of a random surface embedded in a 
three dimensional space, suspended above a planar substrate with which it interacts [IHS]. A detailed specification 
permitted the exact mapping of configurations of the three-dimensional surface onto those of the planar Ising model 
in its Peierls-contour representation. When treated by equilibrium statistical mechanics, this model displayed a phase 
transition; at low enough temperatures, typical configurations indicate that the surface is bound to the substrate. On 
raising the temperature sufficiently, the interface unbinds from the substrate and a layer of the bulk phase intercalates 
between the random surface and the substrate. Although it might be thought, as the authors originally did, that this 
model affords an example of wetting in 2-d [I] , it turns out that the exact critical indices obtained are quite unlike 
those normally associated with wetting. Rather, this model is much more accurately considered as an example of 
either Stransky-Krastanov (SK) [5], or Volmer- Weber (VW) [B] behavior, 
zi. We now go back sixty years to the seminal work of Burton, Cabrera and Frank [7] , who were the first to propose 

that a surface phase transition in a 3-d uniaxial system, say the spin-1/2 Ising model, should display singular behavior 
like the 2-d Ising model. They considered symmetry-breaking boundary conditions for which boundary spins are fixed 
to point up in the upper half space and down in the lower half space. They reasoned that below the 3-d critical 
temperature the upper half space will be in the up magnetized state and the lower half space in the down magnetized 
state. In the mean field approximation, the 2-d lattice of spins at the interface between the oppositely magnetized 
bulk phases is subjected to a mean field that cancels out. Hence this layer should be strictly 2-d in zero magnetic 
field and behave accordingly. This scenario is known to be incorrect; the BCF transition as originally conceived is 
actually a roughening transition [5] with quite different characteristics. The work in [TJ [3] affords examples in which 
the Burton, Cabrera and Frank type of transition is obtained by exact analysis, but in a different scenario and without 
making a mean field approximation. 
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Our first step will be to describe the model as originally formulated [IH3] and to review the results obtained for it. 
We will then indicate why we think the SK or VW scenarios are more appropriate than the wetting one. After that, 
we will point out a number of limitations that were necessary to obtain exact results and then show how Monte Carlo 
simulations can be used to get information when these limitations are relaxed and the exact result route is no longer 
available to us. 

In the spirit of Kossel and Stransky (KS) (9] [10] , configurations of molecules adsorbed on the substrate plane are 
constructed by regarding each molecule as a unit cube the lower side of which meshes with the underlying simple 
square lattice of the substrate, denoted A C Z 2 . Molecular rafts are assembled as close packed arrays of the KS 
cubes. It is clear that, because it has no interior holes, a raft can be described just as well as a simple closed loop on 
A. Loops can meet at vertices of A but not on edges, as this would imply over-counting. The upper faces of the KS 
cubes have height 1. To extend this model further, we permit placing rafts on top of rafts, without overhangs. In this 
way, we encounter loops within loops on A. As an additional restriction, the significance of which will soon become 
apparent, we do not allow any edge of a raft to lie directly above one of any lower raft, thus excluding multiple height 
jumps. This model thus manifests stepped towers erected on the plane; it was termed the "multi-ziggurat" (MZ) [3] 
model in recognition of its architectural precursor in ancient Sumeria. To complete the definition of the MZ model for 
equilibrium statistical mechanics, we must specify the configurational energy. The surface of a collection of molecular 
rafts is composed of two parts, one in contact with the substrate and the other in contact with the surroundings. 
Representing this collection as the equivalent loop form, labelled T, the energy E(T) is given by 

E(T)=rL(T) + (r-e)A 1 (T), (1) 

where r is the surface tension, e is the adhesion energy of a molecule to the substrate, and L(T) is sum of the lengths 
^(7) of the simple closed loops of T, 

L(T) = J2^7) (2) 

The term A\ is the area of contact of the plaquettes with the substrate, which is exactly the same as the "roof area 
because of the ziggurat construction. Thus the second term in Eq. (IT]) is the energetic contribution of plaquettes in T 
with normal perpendicular to the substrate plane. The remainder of the surface tension contribution is given by the 
first term. For stability against detachment, we evidently require that r > — e. We now briefly mention the results 
that can be obtained for this model. Notice that when r = e we have recaptured precisely the planar Ising model, 
so the transition temperature and the free energy singularity are known. This is also the case if e > t; then the first 
layer is completely covered and subsequent rafts are laid on top of this layer as with r = e. In the remaining region 
of stability, characterised by — r < e < t, much less detailed information is available [2]. 

If we are to regard the MZ model as a version of the VW and SK scenarios, then we should point out three 
shortcomings of the model as it stands. The first, which we will not discuss further here, is that we do not allow 
for the energy of mismatch, elastic in origin, between the first raft and the substrate. The second deficiency, the 
examination of which is the main subject of this paper, is the Ehrlich-Schwoebel (ES) |11[ 112] phenomenon. Simply 
stated, multiple height steps should be allowed, but they are energetically disfavored. In the next section, we will 
describe in some detail a model of [13] that is simply related to the MZ model but permits multiple height steps and 
incorporates the ES idea. 

The third problem is the formation of corrals. A corral or hole is a region of lesser height surrounded by a region 
of greater height. The model studied in this paper forbids corrals. Detecting a corral requires non-local information 
and it is unlikely that a physically reasonable equilibrium model would contain long range interactions that forbid 
corrals. On the other hand, dynamical considerations may suppress corrals. The formation of a corral from a plateau 
requires the removal of a molecule that is entirely surrounded by neighbors. Such a molecule will be tightly bound 
and its evaporation suppressed. A corral might also be formed by the sequential deposition of a wall that eventually 
encloses a region. However, surface diffusion would tend to favor more compact structures and suppress the formation 
of walls. Thus, although our model is an equilibrium statistical mechanical model, we believe that the no-corrals rule 
makes it an appropriate model for non-equilibrium surface growth processes. 

From the point of view of equilibrium statistical mechanics, the above height models can be regarded as gen- 
eralizations of the height representation of the 2-d Ising model, which corresponds to Eq. (flj) with r = e. Other 
generalizations are possible. In this work we consider the generalization of [13] that incorporates the ES idea by 
allowing height steps greater than one. Height steps greater than one are given an extra energy proportional to e±. 
We study this multistep height model numerically and find that, along its critical curve, the critical exponents vary 
continuously with t\. Though our model is motivated by surface physics, it appears to be closely related to another 
class of models in statistical physics, the O(n) loop models. In the O(n) loop models the statistical weight depends 



on the total loop length L(T) and also contains a factor n for each simple loop. Thus, the energy E(T) becomes 



E(r)=rL(r) + (lnn)C(r), 



(3) 



where C(T) is the number of simple loops in T and n is referred to as the loop fugacity. On lattices with vertices of 
degree 3, the loops are disjoint and on the honeycomb lattice one has the well-known O(n) loop model introduced 
in [13]. A great deal is known about this O(n) loop model when n < 2. For a given n, there exist three distinct phases: 
a disordered phase with small loops, a densely-packed phase, and a fully-packed phase. The densely-packed and fully- 
packed phases are both critical-i.e., the probability for two points to be on the same loop decays algebraically with 
distance. Between the disordered and the densely-packed phase is a critical curve as a function of n. The singular 
behavior of the O(n) loop model along this critical curve can be described by a set of critical exponents that are 
functions of n. These exponents can be obtained from a mapping to Coulomb-gas theory |15j . Surprisingly, we find 
strong numerical evidence that the multistep height model studied here can also be described by the Coulomb gas 
theory and that there is a one-to-one mapping for universal quantities between n and ei. In the continuum scaling 
limit when the lattice spacing shrinks to zero, critical O(n) loop models and hence presumably also critical multistep 
height models should be conformally invariant and described by Conformal Loop Ensembles — see, e.g., [16] — with 
a parameter value K mapped onto n and t\. 

The plan for this paper is as follows. In Sec.[TTJ we define the multistep height model. In Sec. |III| we list quantities 
of interest for the model. In Sec. IV we describe the numerical methods that are used in our simulations. In Sec. [VJ 
we present our results for the critical behavior of the multistep height model. The paper closes with a discussion in 
Sec.lVl] 



II. THE MULTISTEP HEIGHT MODEL 




FIG. 1: A typical configuration of the multistep height model near the critical temperature for r 
of L=30. The shade of the columns represents their height. 



2, ei = and system size 



The height models proposed in [T3] and studied here generalize the height representation of the Ising model. 
Consider the two-dimensional Ising model on a square lattice in the spin representation with uniform plus-spin fixed 
boundary conditions. There is a one-to-one mapping from spins on the direct lattice to heights on the dual lattice. 
The height of any dual lattice site is the least number of Peierls contours that must be crossed on a path from the 
boundary. This definition leads to several important consequences. First, the magnetization of the Ising model is 



simply the number of even height sites minus the number of odd height sites. Second, there is a constraint that no 
holes or corrals are allowed in the height representation. That is, from every dual lattice site, there exists a path to 
the boundary that is non-increasing in height. Finally, the height representation of the Ising model will only have 
height steps of +1 or — 1. 

To generalize the model, we break the constraint of single height steps by allowing larger height steps at an extra 
energy cost, while still disallowing holes or corrals. The probability for a loop configuration is obtained from the 
energy, E(T), 

E(T)=tL(T) + 6^(1,1). (4) 

Here, T is a collection of simple closed loops defined on the dual lattice. In the Ising model, these form the Pcicrls 
contours, while in the multistep model they can have a more complicated structure. In particular, this model allows 
Peierls contours to overlap, corresponding to larger height steps. L(T), previously defined in Eq. (T2J) , is the sum of 
the lengths of all of the loops (sum of all the edge weights), including now the lengths of all of the overlapping loops. 
In the loop representation r is the energy associated with adding a single edge to T while N(l,l) is defined as L(T) 
minus the length of all edges with weight one (step size one) . ei is the energy cost associated with these larger height 
steps. L(T) adds factors of r linearly in the size of the height step, while N(l, 1) adds factors of e\ linearly in the size 
of the step but only for height steps larger than one. Figure 1 shows a typical configuration of the multistep height 
model near the critical temperature for r = 2, e = 0, and system size L — 30. Note that we have not included the 
term (r — c)Ai(T) in Eq. (fil) in the multistep height. It would be interesting to consider its effect in future studies. 

Equation Q completely specifies the loop representation of the model. The height representation is uniquely 
obtained from the loop representation. The height associated with a lattice site is obtained using the rule that the 
height change across a dual edge is equal to the weight of the bond. The direction of the height change is determined 
by the no-corrals rule. Figure [2] shows a possible edge configuration and the associated heights. In this configuration 
L(T) = 62, N(l, 1) = 6 and E(T) = 62r + 6ei. 

In this paper we simulate multistep height models for several values of e± in the range —1 < ei < oo and r = 2. 
Equation Q with r = 2 and e\ = oo corresponds to the height representation of the Ising model with interaction 
energy J = 1. 

III. MEASURED QUANTITIES 

Order parameters for height models |17H19j can be constructed in analogy to the magnetization for spin models. 
In the height representation of the Ising model, plus spins correspond to even heights while minus spins correspond 
to odd heights. One can define magnetization- like quantities, M n for positive integers n as is done in [19], 

M„ = l£exp(^), (5) 

TV *-^ n + 1 

3 

where h(j) is the height at site j and N is the number of sites. This family of order parameters is motivated by 
the idea that in the rough (high temperature) phase, any height is nearly equally likely to occur so M n — > 0. In the 
smooth (low temperature) phase, a single height will dominate and M n — > 1 as T — > 0. The magnetization in for the 
Ising model in the spin representation is equal to Mi . In the critical region of the multistep model, we find that only 
Mi is useful for the small systems studied here. Heights significantly greater than two do not occur often, therefore 
M n will have strong finite-size corrections for n > 2. Henceforth we use m to represent either the magnetization for 
the spin representation of the Ising model or M\ for height representations. 
We measured the following quantities obtained from the height field h(j): 

Binder cumulant The Binder cumulant U is defined as 

Crossings of the Binder cumulants as a function of temperature for various system sizes are used to locate the 
critical point. 

Order parameter From the finite-size scaling of the order parameter (m) at the critical temperature we obtain the 
exponent ratio /3/v from (m) ~ L~PI V where L is the system length. 



Susceptibility The susceptibility is defined as 

X = /3N(m 2 - (m) 2 ). (7) 

From the finite-size scaling of the susceptibility \ a t the critical temperature we obtain the exponent ratio j/i> 
from x ~ L~ ( >". At the critical temperature for fixed boundary conditions the finite-size scaling relation for x 
is better behaved if the mean is set to zero, the infinite system critical value. Our finite-size scaling fits for 7/1/ 
are made with the subtraction term (m) 2 omitted. 

Specific heat The specific heat c is defined as 

c=§«£ 2 }-<£) 2 )- (8) 

From the finite-size scaling of the specific heat c at the critical temperature we obtain the exponent ratio a/u 
from c - L a / V . 

Bare substrate areal fraction The areal fraction of the bare substrate <f> is defined as 

1 N 
<t>=Nz2 S hU),o- ( 9 ) 

For T <T C there are almost no ad-atoms and « 1 while for T 3> T c , ^«0. The bare substrate is also known 
as the "gasket" in [TS] where its (mean) fractal dimension at criticality is calculated rigorously for conformal 
loop ensembles and agrees with earlier non-rigorous results of Duplantier |20j for the area of connected regions in 
critical O(n) loop models. Following the notation of Ref. [ST] the fractal dimension of these domains is 2 — 0' Jv 
so that the finite-size scaling of <fi is given by 4> ~ L^ 13 l v , 

Height We measured the height at the origin h(0) and the average height of the lattice h = (1/N) J^i h(i). These 
quantities arc expected to diverge logarthmically in the system size. 

IV. NUMERICAL METHODS 

We developed two Monte Carlo algorithms to simulate these height models. For the Ising case, we simulate the spin 
representation using the Wolff algorithm adapted to fixed boundary conditions. The spins are then mapped to heights 
in order to measure properties of the height model. For the multistep height model, there is no spin representation and 
we use a variant of the worm algorithm to sample the collection of closed loops T. The non-locality of the no-corrals 
rule creates significant computational difficulties in implementing the worm algorithm for the multistep height model. 

A. Wolff Algorithm 

The standard Wolff algorithm [22] for the Ising model must be modified for fixed boundary conditions. Our approach 
is to reject clusters that add a boundary spin to the cluster. Here are the steps of the modified Wolff algorithm: 

1. Choose a random site i to initiate the cluster. 

2. Using a breadth-first search, consider all neighbors j of every site in the cluster and propose to add j to the 
cluster. 

• If the neighbor j has the same spin as the cluster add this site to the cluster with probability p = 1 — e~ 2/3 . 

• If the site that has been added to the cluster is part of the fixed boundary, reject the entire cluster. (The 
rejection of clusters that connect to the boundary is the only new feature required for fixed boundary 
conditions. It reduces the efficiency of the algorithm compared to free or periodic boundary conditions.) 

• Continue adding sites to the cluster until either the cluster is rejected because it touches the boundary or 
the growth of the cluster terminates. 

3. Flip the cluster with probability 1. 



4. Collect statistics. 

In order to collect statistics on the height field we must map the spin configuration to a height configuration. We 
do this using a series of breadth-first searches, starting from the boundary. 

1. All the boundary sites are put in the current queue and assigned height zero. All other sites have no assigned 
height. 

2. Do a breadth first search starting from the current queue. Sites are removed from the current queue one at a 
time and all neighbors of the site are tested. Neighbors with the same spin as the current queue and no assigned 
height are added to the current queue and given the same height as the current queue. Neighbors with opposite 
spin and no assigned height value are assigned a height value one greater than that of the current queue and 
added to the future queue. This process is repeated until the current queue is empty. 

3. The sites in the future queue are moved to the current queue and the future queue is emptied. 

4. Steps 2 and 3 are performed iteratively until the heights of the entire lattice have been determined. 

Although the fixed boundary conditions and the spin-to-height mapping slow down the computation, the Wolff 
algorithm is still very efficient for the height representation of the Ising model with fixed boundary conditions. It is 
significantly faster than either the single-spin Metropolis algorithm or the worm algorithm described below for the 
multistep height model. However, the standard worm algorithm for the Ising model 23 might be more efficient for 
fixed boundary conditions but was not used for this study. 

B. Worm Algorithm 

The worm algorithm was developed as a local update algorithm that, nonetheless, almost entirely eliminates critical 
slowing down [23-25]. An additional benefit of using the worm algorithm for the multistep height model is the ability 
to efficiently detect and reject moves that would create corrals. The worm algorithm generates a biased random walk 
that creates directed edges on the dual lattice and samples the closed loop structures T with probabilities associated 
with the energy defined in Eq. Q. 

Height models may be represented either by oriented or unoriented loops. The energy is the same in either case 
but by using oriented loops it is possible to locally encode the direction of height steps and it is more straightforward 
to compute heights and detect corrals. In order to easily satisfy the no corrals rule, the worm algorithm actually 
samples configurations of oriented loops, V . We define height changes across an edge relative to the direction that an 
edge is approached: crossing a right-pointing edge corresponds to an increase in height while crossing a left-pointing 
edge is a decrease in height. (If e is the direction of an edge and v the direction the edge is crossed then the direction 
of the height step is e x i.) As a result counter-clockwise loops enclose mounds while clockwise loops enclose holes 
or corrals and clockwise loops are forbidden by the no-corrals rule. Figure [2] shows an example of an allowed set of 
oriented loops and the corresponding heights. 

The worm algorithm generates properly weighted, oriented, counter-clockwise loop configurations as follows: 

1. Initiate the worm at a random location on the dual lattice. Set the head and tail of the worm at the same 
position on the dual lattice. 

2. Grow the worm. With equal probability, propose to move the head or tail a single step in a random direction 
along the dual lattice. The movement of the head creates a directed edge pointing from the original position of 
the head to its new position. The movement of the tail creates a directed edge pointing from the new position 
of the tail to its original position. Thus, as the head and tail move, the path connecting the head and tail is 
created or destroyed. 

3. For the proposed move to be accepted, two conditions must be met: the move must be accepted on energetic 
grounds and a corral must not be created. 

(a) Calculate the change in energy AE associated with the move and the associated Boltzmann factor e~@ AE . 
An increase in the edge weight from zero to one is provisionally accepted with probability e~^ T while an 
increase from a nonzero value to one higher is provisionally accepted with probability e - "( r+£l ). Decreasing 
the edge weight is always provisionally accepted. These acceptance probabilities insure that the probability 
distribution for loop configurations satisfies detailed balance with respect to the set of worm moves. 



(b) Verify that no corral is created. As noted previously, this means that no clockwise loop is formed. If the 
proposed edge joins two existing edges, then a depth-first search that follows the most clockwise edges is 
carried out. To detect a clockwise loop, we assign every vertex i of the dual lattice an angle 0j. This angle 
is measured from the point at which the search starts. Right (left) turns yield changes in angle of —1 (+1). 
If the search returns to the initial point and has A9i = —4, then it is evident that a clockwise loop has 
been formed and the move is rejected. The search ends by returning to the initial point or when all edges 
connected to the initial point have been explored. 

4. If the edge is provisionally accepted and does not create a corral, the proposed move is accepted and the edge 
is added to V . 

5. Repeat steps 2-4 until the head and tail of the worm meet. It is only when the head and tail meet that the 
edge set is a collection of loops and thus a physical configuration. 

6. Heights are measured when the worm closes and r" is physical. The height of each site is found by scanning 
through the lattice, row by row, and examining the changes in height by counting the edges that are crossed. 
The direction of the edge that is crossed determines whether the height change is positive or negative, while 
the weight of the edge determines the change in height in crossing the edge. All other observables are obtained 
from the heights. 

We find that of these two algorithms, the Wolff algorithm applied to the spin representation followed by a mapping 
to the height representation is the most efficient way to study the height representation of the Ising model. The 
worm algorithm and the Wolff algorithm are comparably efficient as measured in Monte Carlo sweeps, but the worm 
algorithm has a substantial computational overhead associated with checking the orientation of loops. However, the 
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FIG. 2: A configuration of the multistep height model in the oriented loop represenation. Each edge corresponds to a height 
change given by its weight. A double arrow correspond to an edge with weight two, while a single arrow corresponds to an 
edge with weight one. Crossing an edge with an arrow pointing to the right of the crossing direction is a positive height step. 
Numbers indicate the heights of the regions. For this configuration L(T) = 62 and N(l, 1) = 6 and E(T) = 62r + 6ei. 



Wolff algorithm is applicable only in the Ising case where a spin representation is available. The multistep height 
model has no known spin representation and thus it is necessary to work either in the height representation or the loop 
representation. Both the worm and Wolff algorithms are far more efficient than a single site height update algorithm 
such as the single-spin Metropolis algorithm or loop algorithms that update loops locally. 

For the worm algorithm time can be measured in the number of worms attempted. A worm attempt consists 
of choosing a random site and growing the worm from this site until it either closes to form a loop or disappears. 
Measured in units of worm attempts, we found that observables approach their equilibrium values exponentially with 
a time scale that is approximately independent of both the observable and the value of e\ . This equilibration time 
scale grows with system size but, for the largest system L = 199 it is less than 7500. For all system sizes we discard the 
first 10 5 worm attempts to achieve equilibration before collecting data. Data is collected every 200 worm attempts. 

Errors in the measurement of observables are estimated using the blocking method |26] . Each simulations consist 
of 200,000 data collection sweeps (4 x 10 7 worm attempts) after the initial equilibration. The data is divided into 
200 blocks each containing 1,000 measurements and errors are obtained from the standard deviation between blocks. 
The critical temperature is estimated by extrapolating the crossings of the Binder cumulants to infinite system size as 
described below. The error in the critical temperature was obtained by considering power law fits of critical observables 
near this extrapolated critical temperature. The range of temperatures for which the power law fit was a good fit 
determined the error in the critical temperature. Errors in extrapolated quantities such as the critical exponents 
are estimated from the range of fit parameters found by varying the critical temperature over its uncertainty range. 
Statistical errors in the fits needed to obtain exponents are much smaller, typically 10% to at most 25% of the 
uncertainty resulting from the uncertainty of the critical temperature. 

V. RESULTS 

A. Critical Ising model in the height representation (ei = oo) 

To understand critical properties of the height representation of the Ising model, the Wolff algorithm with fixed 
boundary conditions, as described in section [IV A[ is used along with finite-size scaling. We simulated the system at 
the critical temperature T c = 2/ln(l + y/2) for 12 sizes ranging from L = 19 to L — 199 to determine the scaling 
behavior of the average height per site. Odd system sizes are used so that the domain has a unique center. 
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FIG. 3: The average height h versus L for the Ising model. The solid line is the best fit to the form given in Eq. ( 10 \ where 
the prefactor to the logarithm B = .0467. 

Figure [3] shows how the average height per site scales with system size. The fitting function is 

h = A + Bhx(L)(l + T ). (10) 
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The 1/L correction to scaling is suggested by the fixed boundary conditions. It can be difficult to distinguish between 
logarithms from small power laws, but the data clearly shows that this is a logarith mic function. A three parameter 
power law fit for h yields a Q-value less than .0001, while the logarithmic fit of Eq. |l0| ) has a Q value of 0.34. The 
prefactor of the logarithm is B — 0.0467 ± .0015. This prefactor is the same within error bars for h as it is for h(0). 

Using finite-size scaling we also obtain values for (3/v, /3' /v and y/v as shown in Table IT] These and other critical 
exponents are obtained from a three parameter fit of the form y = aL b (l + c/L) where, again the linear correction to 
scaling is suggested by the fixed boundary conditions and yields reasonable fits. Although fixed boundary conditions 
results in larger finite-size corrections than free or periodic boundary conditions we are still able to obtain accurate 
exponent values. 



B. Critical Behavior of the Generalized Height Model (ei < oo) 



The critical temperature of the multistep model depends on t\ . Critical temperatures are obtained from crossings 
of the Binder cumulant defined in Eq. (pi). Figure R] shows the Binder cumulant for various system sizes. Crossings 
of these curves for successive system sizes gives a finite-size estimate of the critical temperature. By plotting the 
crossings vs. 1/L and extrapolating to 1/L —} 0, we find an estimate of the critical temperature in the thermodynamic 
limit, as seen in Fig. [5] This method is used for various values of ei. The results are shown in Table IT] 
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FIG. 4: (color online)The Binder cumulant U is shown as a function of temperature T for various system sizes and ei = 0. 
System sizes increase from top to bottom on the right side of the graph and from bottom to top on the left side. Finite-size 
estimates of the critical temperature are obtained from the crossings of the curves for successive system sizes. 



Using these values for the critical temperature, simulations are carried out at and near the measured critical 
temperature for 12 system sizes from L = 19 to L = 199 to determine the finite-size scaling behavior of each of the 
quantities of interest. 

We find that the multistep model has critical behavior that is dependent on the value of €\. We study the exponents 
(3/v, 7/V, a/v, and (3' /v, as well as the logarithmic prefactor B of the average height per site. Table [TJ lists the 
resulting fits for each value of t\. The errors in the specific heat measurements are too large to discriminate between 
a logarithmic scaling law (a/v = 0(log)) and a small power law a/v > 0. Both fits are plausible and yield similar 
goodness of fit values. The values of a/v in Table [TJ assume a power law fit. 

It is possible that the observed continuous variation in the exponents shown in Table |l] reflects finite-size corrections 
and that the number of universality classes is small, perhaps one or two. However, we believe that the evidence points 
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FIG. 5: The crossings of U are plotted vs. 1/L for systems between L 
the smaller of the system sizes. 



60 and L = 160. The L value used for each point is 



ei 


T c 


P/v 


■y/v 


a/v 


l3'/v 


B 


-1 


2.0875(10) 


.065(3) 


1.868(8) 


.33(3) 


.0317(15) 


.031(2) 





2.1685(10) 


.075(2) 


1.85(1) 


.26(2) 


.034(1) 


.0326(15) 


1 


2.2065(10) 


.083(2) 


1.83(1) 


.21(2) 


.038(2) 


.038(2) 


2 


2.230(1) 


.0925(20) 


1.815(5) 


.17(2) 


.042(1) 


.040(3) 


oo (simulation) 
oo (exact) 


2.2686(8) 
2.269185. . . 


.127(4) 
.125 


1.75(2) 
1.75 





.053(2) 
.05208... 


.0467(15) 
.04594... 



TABLE I: Critical temperature T c , critical exponents j3/v , y/v, a/v and /3' /v, and the prefactor of the logarithmic scaling of 
the average height B for each simulated value of ex- The Ising model corresponds to ei = oo. Ising model results were obtained 
from the Wolff algorithm. Exact results for the Ising model are presented for comparison. 



to the conclusion that the multistep height model describes a one parameter family of universality classes parameterized 
by ei. In two dimensions, there are several parameterized statistical mechanical models that have continuously varying 
critical exponents. These include the g-state random-cluster model in terms of parameter q £ [0,4] [57], the O(n) 
loop model for n £ [—2, 2], the Ashkin- Teller model in the associated coupling strength, as well as the 6- or 8- vertex 
models. In all of these systems, the critical exponents can be expressed in terms of a single parameter-the coupling 
constant g in the Coulomb-gas model. Thus, it is plausible that g is also sufficient to describe the critical exponents 
of the multistep height model. 

We note that the multistep height model and the O(n) loop model [21 [3TJ [551 US] are similar in several ways. Both 
models are formulated in terms of loops and contain the Ising model as a special case. Both models have statistical 
weights that depend in the same way on the total loop length, L(T). The O(n) loop model is typically denned on 
a honeycomb lattice and requires Eulerian loops so that loops do not overlap and the statistical weight has a term 
n c ( r ) where C(T) is the number of simple loops. When n = 1 the model reduces to the Ising model. In Ref. |28j . 
by making use of the fact that the critical O(n) loop model is equivalent to the tricritical q = n 2 Potts model, the 
critical exponents for the magnetization, the Ising-spin domain, the susceptibility, and the specific heat were related 
by simple formulae to the Coulomb gas coupling. In terms of the exponents f3/v,/3' /i/,j/i/, and a/v, these formulas 
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read [JJ3 CUJ MM] ■ 



P/v = 


i 
9 




fflv = 


(g-2)(6- 


^5) 


8.9 




l/v = 


4.9 - 12 
9 




ajv = 


r 32 
6 , 





(11) 

(12) 
(13) 
(14) 



with g <E [4,6]. Equation (13) is related to Eq. (11) by j/p = 2 — 2f3/v. Later in this section we will discuss the 
relation between the constant B in Equation ( 10 1 and the Coulomb gas parameter g. 

We use our measurements of /3/v, P'/is, "f/v, and ajv to estimate g for each value of e\. Estimates of g from 
these exponents agree well with each other and a combined estimate can be obtained from a weighted average of the 
estimates from each of the four critical exponents. The weight assigned to exponent i is l/Sf where 5i is the error 
in the estimate of the corresponding exponent. Table ITT] lists the value of g, calculated by performing a weighted 
average of each of the four measured exponents together with the goodness of fit Q for combining the values from each 
exponent into a single value of g. The goodness of fit is calculated assuming the critical exponents are independent 
and normally distributed with a standard deviation given by Si. These assumptions are not quantitatively correct so 
the Q values cannot be interpreted as the probability of finding \ 2 larger than the fit value. Nonetheless, the fact 
that the Q values are large suggests that a single g simultaneously predicts all of the exponents within their error 
ranges. The existence of a single g consistent with all the critical exponents supports the idea that the multistep 
height models parameterized by e\ are indeed in the universality classes of the Coulomb-gas model parameterized by 
g. The last column in Table In lists the value of the loop fugacity n in the O(n) loop model for the value of g, which 
is obtained from the formula [151 12"5] 



n = — 2cos(7rg/4). 



(15) 



The relation between n and ei for the four finite values of e\ that were simulated is close to linear and extrapolates 
to n = 1/2 for t\ = —2, which is the limit of stability of the multistep height model. 



ei 


ff(flt) 


Q 


ff(conj) 


n(fit) 


-1 


5.629(10) 


.759 


5.6261 


.574(15) 





5.579(8) 


.986 


5.5759 


.649(11) 


1 


5.535(9) 


.802 


5.5297 


.714(13) 


2 


5.488(7) 


.731 


5.4890 


.783(10) 


oo (measured) 


5.324(16) 


.977 


5.3333 


1.013(43) 


oo (exact) 


16/3 


- 


- 


1 



TABLE II: For each simulated value of ei, the best fit values of the Coulomb gas coupling g(fit), the goodness of this fit Q, the 
conjectured value of the Coulomb gas coupling g(conj) obtained from Eqs. ( 15 1 and ( 16 1, and the loop fugacity Ti(fit) obtained 
from (?(fit) and Eq. ( 15 I. 



We find that the simple relation 



tanh[(ei +2)/2tt] + 1 



(16) 



is a good fit to the data in Table III] The function is a guess based on the requirement that it yields n = 1/2 for 
ei = —2 and approaches n — 1 exponentially as ei — > oo as might be expected since t\ exponentially suppresses loop 
weights greater than one. Figure [6] shows the data in Table [IT] together with the conjectured relation, Eq. (16). It 
would be interesting to simulate larger values of t\ to test this conjecture beyond the linear regime. If the conjecture 
holds, it suggests the possibility of an exact mapping from the multistep height model to a known 2-d model. 

There are also exact results for the Coulomb gas model for the universal prefactor of the logarithmic growth of the 
height at the origin. It is straightforward to show that the same prefactor holds for bo th t he height at the origin and 
the average height, h. The prediction for the logarthmic prefactor B defined in Eq. (10) from [32], which has been 
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FIG. 6: The loop fugacity n of the O(n) loop model vs. e\. The curve is obtained from Eq. ( 16 1 and the data is from Table |TT| 



proved for conformal loop ensembles |16j . is 



B= ^ ^COt(W4). 

■Kg 



(17) 



If we use the values of g given in Table|lT|to compute B using Eq. ( fl7| ) and compare with the measured values in Table 
III we find reasonable agreement. The measured values are larger than the predicted values by about twice the quoted 
error except for the Ising case where the two values differ by about the quoted error. Since the measurement of B is 



likely to have substantial systematic errors that are not included in the 1/L correction in Eq. (10 1, we believe that 
the results for B are consistent with the Coulomb gas predictions and provide additional evidence that the multistcp 
height models are in the universality classes of the Coulomb gas/conformal loop ensemble. 



VI. DISCUSSION 



We have analyzed a height model that generalizes the height representation of the Ising model by allowing height 
steps greater than unity with an energy cost parameterized by t\. Our data supports the hypothesis that the critical 
exponents and prefactor B depend continuously on e%. Although we believe that the system has continuously varying 
critical exponents, it is conceivable that there are only one or two universality classes and that the ei dependent 
exponents reflect a finite-size crossover. The hypothesis of continuously varying exponents is strengthened by the 
relationship between the measured exponents that allows a mapping via the Coulomb gas parameter g to the O(n) 
loop model for n in the range 1/2 < n < 1. We conjecture an exact relationship between e\ and g or n that merits 
further investigation. 

In the multistcp height model the energy of a height step is a linear function of its magnitude. It would be interesting 
to study other energy functions. For example, the energy of a height step could be quadratic in its magnitude. We 
suspect that this would also yield models in the Coulomb gas universality classes. 

A central feature of the model studied here is the no-corrals rule. It would be interesting to study a model with the 
same energetics but without forbidding corrals. Height models that allow corrals are generally referred to as "solid-on- 
solid" models. If the energy for height steps greater than one is quadratic and corrals are allowed the corresponding 
discrete Gaussian model is known to be in the 0(2) universality class [33] ■ The body-centered solid-on-solid (BCSOS) 
model restricts height steps to one and has an energy that is a sum of next-nearest neighbor height differences. The 
BCSOS model is exactly solvable and is also in the 0(2) universality class [34j|35j. In the loop representation, allowing 
corrals corresponds to oriented loops with both orientations allowed. In the case of oriented loops on the honeycomb 
lattice with ej = oo, the loops are disjoint and the orientation degrees of freedom can be simply integrated out. This 
yields a statistical weight of 2 for each loop, and hence maps to the 0{n) loop model with n = 2. On this basis, it 
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is tempting to speculate that for < e\ < oo, these oriented loop models will map to 0(n) loop models with loop 
fugacity in the range 1 < n < 2. 

Acknowledgments 

We acknowledge useful discussions with Nikolay Prokofiev. J. M. and M. D. were supported by NSF DMR-0907235. 
C.N. was supported by NSF OISE-0730136 and NSF DMS-1007524.. Y. D. was supported by the NSFC under Grant 
No. 10975127. 



[1 
[2 
[3 
[4; 

[5 
[6 
[7 
[8 
[9 
[10 

[11 
[12 

[13 

[14 
[15 
[16 
[17; 
[18 
[19 
[20 
[21 
[22 
[23 
[24 

[25 
[26 
[27 
[28 
[29 
[30 
[31 
[32 
[33 

[34 
[35 



D. B. Abraham and C. M. Newman, Phys. Rev. Lett. 61, 1969 (1988). 

D. B. Abraham and C. M. Newman, Journal of Statistical Physics 63, 1097 (1991). 

D. B. Abraham, L. Fontes, C. M. Newman, and M. S. T. Piza, Phys. Rev. E 52, R1257 (1995). 

D. B. Abraham, in Phase Transitions and Critical Phenomena, edited by C. Domb and J. L. Lebowitz (Academic, London, 

1986), vol. 10. 

I. N. Stranski and L. Krastanow, Akad. Wiss. Wien Math.-Natur. Kl. lib 146, 797 (1939). 

M. Volmer and A. Weber, Z. Phys. Chem. 119, 277 (1926). 

W. K. Burton, N. Cabrera, and F. C. Frank, Philos. Trans. Roy. Soc. London A243, 299 (1951). 

J. D. Weeks, G. H. Gilmer, and H. J. Leamy, Phys. Rev. Lett. 31, 549 (1973). 

W. Kossel, Nachr. Ges. Wiss. Gottingen Math. Phys. KLS 135, 135 (1927). 

I. N. Stranski, Z. Phys. Chem. 136, 259 (1928). 

G. Ehrlich and F. C. Hudda, J. Chem. Phys. 44, 1039 (1966). 

R. L. Schwoebel, J. Appl. Phys. 40, 614 (1969). 

D. B. Abraham and C. M. Newman, EPL (Europhysics Letters) 86, 16002 (2009). 

E. Domarry, D. Mukamel, B. Nienhuis, and A. Schwimmer, Nuclear Physics B 190, 279 (1981). 
B. Nienhuis, J. Stat. Phys. 34, 731 (1984). 

O. Schramm, S. Sheffield, and D. Wilson, Communications in Mathematical Physics 288, 43 (2009). 

J. R. G. de Mendonca, Phys. Rev. E 60, 1329 (1999). 

G. Mazzeo, G. Jug, A. C. Levi, and E. Tosatti, Journal of Physics A: Mathematical and General 25, 967 (1992). 

U. Alon, M. R. Evans, H. Hinrichsen, and D. Mukamel, Phys. Rev. E 57, 4997 (1998). 

B. Duplantier, Phys. Rev. Lett. 64, 493 (1990). 

Q. Liu, Y. Deng, and T. M. Garoni, Nuclear Physics B 846, 283 (2011). 

U. Wolff, Phys. Rev. Lett. 62, 361 (1989). 

N. Prokof'ev and B. Svistunov, Phys. Rev. Lett. 87, 160601 (2001). 

N. Prokof'ev and B. Svistunov, in Understanding Quantum Phase Transitions, edited by L. D. Carr (Taylor & Francis, 

Boca Raton, 2010). 

Y. Deng, T. M. Garoni, and A. D. Sokal, Phys. Rev. Lett. 99, 110601 (2007). 

M. E. J. Newman and G. T. Barkema, Monte Carlo Methods in Statistical Physics (Oxford University Press, 2007). 

C. M. Fortuin and P. M. Kasteleyn, Physica 57, 536 (1972). 
B. Nienhuis, Phys. Rev. Lett. 49, 1062 (1982). 

Y. Deng, T. M. Garoni, W. Guo, H. W. J. Blote, and A. D. Sokal, Phys. Rev. Lett. 98, 120601 (2007). 

H. Saleur and B. Duplantier, Phys. Rev. Lett. 58, 2325 (1987). 

B. Duplantier and H. Saleur, Phys. Rev. Lett. 63, 2536 (1989). 

J. Cardy and R. Ziff, Journal of Statistical Physics 110, 1 (2003). 

B. Nienhuis, in Phase Transitions and Critical Phenomena, edited by C. Domb and J. L. Lebowitz (Academic Press, 

London, 1987), vol. 11, p. 1. 

H. van Beijeren, Phys. Rev. Lett. 38, 993 (1977). 

R. J. Baxter, Exactly Solved Models in Statistical Mechanics (Dover Publications, 2008). 



